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ABSTRACT 

The importance of the radiative feedback from massive black holes at the 
\ centers of elliptical galaxies is not in doubt, given the well established relations 

among electromagnetic output, black hole mass and galaxy optical luminosity. 
We show how this AGN radiative output affects the hot ISM of an isolated 



O 



elliptical galaxy with the aid of a high-resolution hydrodynamical code, where 
the cooling and heating functions include photoionization plus Compton heating. 
We find that radiative heating is a key factor in the self-regulated coevolution 
of massive black holes and their host galaxies and that 1) the mass accumulated 
by the central black hole is limited by feedback to the range observed today, 
and 2) relaxation instabilities occur so that duty cycles are small enough (<0.03) 
to account for the very small fraction of massive ellipticals observed to be in 
the "on" -QSO- phase, when the accretion luminosity approaches the Eddington 
luminosity. The duty cycle of the hot bubbles inflated at the galaxy center during 
major accretion episodes is of the order of ^0.1 — 0.4. Major accretion episodes 
caused by cooling flows in the recycled gas produced by normal stellar evolution 
trigger nuclear starbursts coincident with AGN flaring. During such episodes the 
central sources are often obscured; but overall, in the bursting phase (1<^<3), 
the duty-cycle of the black hole in its "on" phase is of the order of percents and it 
is unobscured approximately one-third of the time. Roughly half of the recycled 
gas from dying stars is ejected as galactic winds, half is consumed in central 
starbursts and less than 1% is accreted onto the central black hole. Mechanical 
energy output from non-relativistic gas winds integrates to 2.3 x 10 59 erg, with 
most of it caused by broadline AGN outflows. We predict the typical properties of 
the very metal rich poststarburst central regions, and we show that the resulting 
surface density profiles are well described by Sersic profiles. 
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Subject headings: accretion, accretion disks — black hole physics — galaxies: 
active — galaxies: nuclei — quasars: general — galaxies: starburst 

1. Introduction 

Supermassive black holes (SMBHs) at the centers of bulges and elliptical galaxies have 
played an important role in the processes of galaxy formation and evolution (e.g., see Silk & 
Rees 1998; Fabian 1999; Burkert & Silk 2001; Cavaliere & Vittorini 2002; King 2003; Wyithe 
& Loeb 2003; Haiman, Ciotti & Ostriker 2004; Granato et al. 2004; Sazonov et al. 2005; 
Murray, Quataert & Thompson 2005; Di Matteo, Springel & Hernquist 2005; Begelman & 
Nath 2005; Hopkins et al. 2006; Croton et al. 2006), as strongly supported by the remarkable 
correlations found between host galaxy properties and the masses of their SMBHs (e.g., see 
Magorrian 1998, Ferrarese & Merritt 2000, Gebhardt et al. 2000, Yu & Tremaine 2002, 
McLure & Dunlop 2002, Graham et al. 2003). 

An additional very severe issue that AGN feedback likely addresses is that of the "cooling 
flow problem" (e.g., Xu et al. 2002, Peterson & Fabian 2006), which in elliptical galaxies 
is at least as serious as the heavily analyzed cooling flow problem in clusters (because the 
radiative cooling times are an order of magnitude shorter, ~ 10 7 - 5 yr vs. ~ 10 9 yr). Also, 
a rarely discussed problem is that the recycled gas (primarily from red giant winds and 
planetary nebulae) of the evolving stellar population expected to accumulate within galaxies 
over cosmic time, and responsible for the cooling flow, is orders of magnitude larger than 
the mass of either the central SMBH or the resident diffuse gas. Quantitatively, the mass 
return rate from evolving stars in elliptical galaxies can be estimated as 

M*(t) ~ 1.5 x — jjy t^ 3 M yr-\ (1) 
iu -Lb© 

where L B is the present galaxy blue luminosity and t 15 is time in 15 Gyr units (Ciotti et al. 
1991, hereafter CDPR, see also Sect. 2.2). Thus, the total mass return would accumulate 
onto the central SMBH a mass by far too large compared with the observed SMBH masses 
if a long lived cooling flow occurred. Young stellar populations observed in the body of 
ellipticals also cannot account for the total mass released, and alternative forms of cold mass 
disposal (such as distributed mass drop-out/star formation), are not viable solutions (e.g., 
Binney 2001). In addition to this mass disposal problem, also the X-ray luminosity Lx of 
low-redshift elliptical galaxies is inconsistent with the standard cooling flow model. In fact, 
low-redshift elliptical galaxies with optical luminosity Lb~3x 1O 1O L show a significant range 
in the ratio of gas-to-total mass at fixed L B , with tabulated values ranging from virtually zero 
up to few % (e.g., Roberts et al. 1991), and most of that is seen in X-rays with temperatures 
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close to the virial temperatures of the systems (~ 10 6 7 K, e.g. O'Sullivan, Ponman & Collins 
2003). 

A (partial) solution to these problems was proposed by D'Ercole et al. (1989) and 
CDPR, by considering the effect of SNIa heating of the galactic gas, and exploring the time 
evolution of gas flows by using hydrodynamical numerical simulations. Subsequent, more 
realistic galaxy models (with updated rates of SNIa derived by direct counts from optical 
observations) were explored by Pellegrini & Ciotti (1998). It was found that while SNIa input 
sufficed for low and medium-luminosity elliptical galaxies to produce fast galactic winds, the 
inner parts of more massive spheroids would nevertheless host inflow solutions similar to 
cooling flows. This is because, while the number of SNIa per unit optical luminosity is 
expected to be roughly constant in ellipticals, the gas binding energy per unit mass increases 
with galaxy luminosity. However, it is not likely that SNIa's can provide the entire solution 
to this problem, because the distribution of the energy input is not concentrated enough 
to balance the observed gas cooling rate (since this scales as p 2 , the required heating rate 
must be very large in the very central regions). Note also that the SNIa rate is independent 
of the current thermal state of the X-ray emitting plasma, therefore SN heating cannot act 
as a self-regulating mechanism. Finally, as already recognized by CDPR, the mass budget 
problem would still affect medium-large galaxies, putative hosts of luminous cooling flows. 
Thus, a concentrated feedback source is a very promising solution for a variety of problems, 
and the central SMBH is the natural candidate, by its mass and by its location, through a 
combination of mechanical and radiative feedback mechanisms (e.g., see Fabian, Celotti & 
Erlund 2006, and references therein). 

Some calculations have allowed for a physically motivated AGN feedback (e.g., see 
Binney & Tabor 1995; Ciotti & Ostriker 1997, 2001, hereafter C097, CO01; Omma et 
al. 2004; Ostriker & Ciotti 2005; hereafter OC05, Churazov et al. 2005), and the computed 
solutions are characterized by relaxation oscillations (Ostriker et al. 1976; Cowie, Ostriker & 
Stark 1978). Energy output (radiative or mechanical) from the central SMBH pushes matter 
out, the accretion rate drops precipitously and the expanding matter drives shocks into the 
galactic gas. Then the resulting hot bubble ultimately cools radiatively (it is thermally 
unstable) and the resulting infall leads to renewed accretion; the cycle repeats. Among the 
computed models which studied the interaction between AGN feedback and galactic cooling 
flows, those of C097 and CO01 focused on the effects of radiative heating on galactic gas flows. 
In fact, if one allows the radiation emitted from the accreting SMBH to interact with and heat 
the galactic gas, one solves the cooling flow problem in elliptical galaxies, and the feedback 
produces systems that are variable but typically look like normal ellipticals containing hot 
gas. They sometimes look like incipient cooling flows and rarely, but importantly, appear like 
quasars. Interestingly, observations seem to support this scenario (e.g., Russell, Ponman, & 
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Sanderson 2007). 

In C097 and CO01, however, a major uncertainty remained about the typical QSO 
spectrum to adopt, in particular the high energy component of that spectrum, which is most 
important for heating the ambient gas. Thus, a simple broken power law was adopted for 
the spectrum with a range of possible values of the Compton temperature - from 10 7 2 K 
to 10 9 ' 5 K - with most of the emphasis of the paper being on the higher temperatures. 
Subsequent work by Sazonov, Ostriker & Sunyaev (2004, see also Sazonov et al. 2007), which 
assessed the full range of observational data of AGNs and computed their Spectral Energy 
Distribution, concluded that the typical equilibrium radiation temperature was narrowly 
bounded to values near 10 7,3 K, i.e., of the order of 2keV. This value, although it is at the 
lower end of the range adopted by C001, is still above the virial temperature of all galaxies 
and, most importantly, well above the central temperature of the cooling flow gas. As noted 
by Sazonov et al. (2005), there is a rather large compensating effect also not included in 
CO01: gas heated by radiation with a characteristic temperature near 10 7 K is heated far 
more effectively by absorption in the atomic lines of the abundant metal species than by 
the Compton process. In particular, Sazonov et al. (2005) provide a fitting formula for 
the Compton plus photoionization and line heating/cooling that we implemented into our 
numerical code, together with additional physics that was missing in C097 and CO01; a 
very preliminary exploration of the new models can be found in OC05. Consistently with a 
second generation of metal rich stars observed in recent SDSS surveys (e.g., Fukugita et al. 
2004; Nolan, Raychaudhury, & Kaban 2007), we now also allow for star formation, which is 
found of primary importance during the first few Gyr, as it consumes a large fraction of the 
cooling gas in a central starburst (e.g., Reuland et al. 2007). 

In this paper we address the coevolution of a SMBH and of the ISM of the host galaxy, 
and argue that the AGN/starburst feedback effects can be the answer to the triple problems 
of (1) the suppression of cooling flows within galaxies, (2) the large scatter in their hot gas 
mass at fixed optical luminosity, and (3) disposition of the recycled gas. It is found that 
the intermittencies described in C097 and CO01 are confirmed also with the more accurate 
treatment of the radiation field, and these have two consequences: they cause nuclear star- 
bursts and they feed the central SMBH in what we observe as coincident AGN/starburst 
episodes. In addition, since the fuel is in fact proportional to the evolving stellar mass, one 
maintains - on average - a central SMBH mass and a younger stellar population mass that 
are both proportional to the galactic stellar mass, with the energetic feedback from central 
SNII and AGN producing the energy input needed to keep the bulk of the gas in a state of 
low density and high temperature. For simplicity in this paper we just present a summary of 
the main properties of a typical model, leaving to following papers the detailed description 
of specific features of observational and theoretical interest. 
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The paper is organized as follows. In Section 2 we describe how the galaxy models 
adopted in the simulations are built, the details of the input physics, and their numeri- 
cal implementation. In Section 3 the time evolution of a representative model galaxy is 
illustrated in detail. Finally, in Section 4 we discuss the main results obtained. 

2. The models 

The galaxy models and the input physics adopted for the simulations have been sub- 
stantially improved with respect to C097 and CO01: in the following we describe them in 
detail for future reference. 



2.1. Structure and internal dynamics 

In CO01 the galaxy models utilized a King (1972) stellar distribution plus a quasi- 
isothermal dark matter halo, in line with the models then used for cooling-flow studies. 
However, the existence of large cores of constant surface brightness has been clearly ruled 
out (e.g., see Jaffe et al. 1994, Faber et al. 1997, Lauer et al. 2005), as HST observations 
have shown how the central surface brightness profile is described by a power-law as far in as 
it can be observed, i.e. to ~ 10 pc for Virgo ellipticals. Following Pellegrini & Ciotti (1998), 
we then adopt a stellar density distribution described by the more appropriate Hernquist 
(1990) model 

^* 2tt r(r + r*) 3 

Optical (e.g., Saglia et al. 1993; Bertin et al. 1994; Cappellari et al. 2006, Douglas et al. 
2007) and X-ray (e.g. Fabian et al. 1986, Humprey et al. 2006) based studies typically find 
that luminous matter dominates the mass distribution inside the effective radius R c , while 
dark matter begins to be dynamically important at 2 — 3R C , with common values of the total 
dark-to-luminous mass ratio 1Z = M^/M* in the range 1^71^6. Theoretical (e.g., Dubinski & 
Carlberg 1991; Ciotti & Pellegrini 1992; Navarro, Frenk & White 1997; Fukushige & Makino 
1997) and observational (e.g., Treu et al. 2006) arguments support the idea that similarly 
to the stellar profiles, also the radial density distribution of the dark halos is described 
by a cuspy profile with central spatial density increasing as r -1 : consistently, for the dark 
matter halo we also adopt the density distribution in equation (2), where M h = TZM* and 
7"h = (3r* are the halo total mass and scale-length; dynamical and phase-space properties of 
the resulting two-component Hernquist models are given in Ciotti (1996). 

In order to use realistic galaxy models, the mass distribution is determined as follows. 



(2) 
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We fix a value for the projected central velocity dispersion a, and we determine the galaxy 
present day model blue luminosity L B and effective radius R c from the Faber-Jakson 

and the Fundamental Plane 

\ogR c = A\oga + B\ogI c + C (4) 

relations adopted in CDPR, where I e = Lb/(2itRI). The scale-length r* of the stellar 
distribution (2) is then given by r* ~ i? c /1.8153 (Hernquist 1990). 

Having fixed a and r*, we determine the galaxy stellar mass and halo properties such that 
<7 a , the model aperture velocity dispersion within R e /S (obtained by solving and projecting 
the Jeans equations), coincides with a. For globally isotropic two-component Hernquist 
models 

where 

"^=< i+ ^()/S- 1 ) 1 (6) 

and 7?.o.5 is the dark-to-visible mass ratio within the stellar half-mass radius (Pellegrini & 
Ciotti 2006); note that we do not consider the effect of M B h on <r a , in accordance with 
estimates of the SMBH sphere of influence radius (e.g., Riciputi et al. 2005). For chosen 
values of 1Z and 7£ .5 we obtain M* and the stellar mass-to-light ratio Y* = M*/L B from 
equations (5)-(6). In the initial conditions we then expand the scale radius r* thus determined 
by a factor 1.5 (while maintaining fixed the dark matter halo distribution and total mass), 
in order to allow for the subsequent shrinkage of newly formed stars in the galaxy central 
regions (see Sect. 2.3). As discussed in the Conclusions, we do not attempt at this stage 
to properly follow the dark matter halo contraction, nor the modifications of the velocity 
dispersion profile consequent to star formation. 

An important ingredient in the energetics of the gas flows, namely the thermalization 
energy of the stellar mass losses, depends on the radial trend of the stellar velocity dispersion 
(see Sect. 2.7) which, for the isotropic case is given (via the Jeans equations) by 



2 2,,, ^ . GM*M BU 



p*a* = p*<7*(Mbh = 0) 



61n(1 + iV (1 + 2s)(6s2 + 6s - i; 



sj 2s 2 (l + s)^ 



(7) 



where s = r/r*, and the first term at the r.h.s. is the 1-dimensional stellar velocity dispersion 
without the contribution of the SMBH (Ciotti, Lanzoni & Renzini 1996). 
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2.2. Stellar passive evolution: SNIa rate and stellar mass losses 

The stellar mass loss rate and the SNIa rate associated with the initial stellar distribution 
are the main ingredients driving evolution of the models. In the code the stellar mass losses 
- the source of fuel for the activity of the SMBH - follow the detailed prescriptions of the 
stellar evolution theory, while for quick calculations the approximation given in equation (1) 
can be used. Over the whole galaxy 

M* = IMF(M to )|Mto|AM, (8) 

where the initial mass function IMF is a Salpeter law (normalized as described in CDPR), 
and the turn-off mass (in M ) of stars at time t (in Gyrs) is 

logM TO = 0.0558(logt) 2 - 1.338 log t + 7.764. (9) 

Finally 

AM = f M to - M fin (M TO ) = 0.945M TO - 0.503, (M TO < 9M ), 

\ AM = M TO - 1.4M , (M TO > 9M ), [ ' 

(CDPR). Recently, updated formulas have become available (Maraston 2005), but in the 
present context the modifications on the flow evolution are minor (Pellegrini & Ciotti 2006) 
and so for continuity we maintained the treatment of past works. Also, SNIa rates in 
early-type galaxies have been carefully reanalyzed, and current estimates of the rate in the 
local universe agree on 0.32/i 2 SNu (where 1 SNu = 1 SNIa per century per 1O 1o Lb and 
h = iJo/100 km s" 1 Mpc" 1 ; e.g., see Cappellaro, Evans & Turatto 1999; Mannucci et al. 
2005) so that, following CDPR, we parametrize the time evolution of the SNIa rate as 

flsN ( ( ) = 0.32xlO-^ SN ^( I ^)" yr~\ (11) 

where the coefficient $sn allow for different choices in the present-day SNIa. Assuming 
for each supernova event an energy release of = 10 51 erg, a fraction 7^sn of which is 
thermalized in the surrounding ISM, the energy input per unit time over all the galaxy body 
is given by 

L SN (t) = 1.015 x 10 31 /A? SN ?7 SN -^ ( 137 f Gyr ) erg 8 "'; ( 12 ) 

in this paper we restrict to the case -#sn = 1 and h = 0.75. As is well known, the specific 
value of s is a critical ingredient in the model evolution. In fact, when s^l.3 the flow evolves 
from wind to inflow, while the opposite is true for s< 1.3. This was especially true in CDPR, 
C097 and CO01 models. However, the specific value of s is less important in affecting the 
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evolution of gas flows in cuspier galaxy models with a somewhat reduced SNIa present-day 
rate, as those here explored, that preferentially host "partial wind" solutions (Pellegrini & 
Ciotti 1998). Thus, here we restrict to the currently favoured s — 1.1 value (Greggio 2005), 
even though a more complicate time dependence than a simple power-law seems possible 
(e.g., Matteucci et al. 2006; Neill et al. 2007). 

Besides energy, supernovae provide also mass. We assume that each SNIa ejects 1.4M Q 
of material in the ISM, so that the total rate of mass return from the aging initial stellar 
population at each place in the galaxy is 

^ = (a* + « S n)p*, (13) 

where asN(^) = 1-4M Rsu(t)/M* and a*(t) = M*{t)/M* are the specific mass return rates. 
With these definitions, the SNIa kinetic energy injection per unit volume in the ISM can be 
written as 

TP TP "^ SN TP a SN(t)p* n ,, 

Ei = ^sn^sn^-P* = ?7sn£sn 1 AMq ■ ( 14 ) 



2.3. Star formation, SNII heating and starburst properties 

At variance with CO01, in the present study we allow for star formation, which cannot be 
avoided when cool gas accumulates in the central regions of elliptical galaxies. In particular, 
we compute the star formation rate at each radius r from the equation 

Pj = , 7f orm = max(r coo i,r dyn ), (15) 

Tform 

where p is the local gas density, reform = 0.03 0.4 (e.g., see Weinberg, Hernquist, & Katz, 
2002; Cen & Ostriker 2006), and the associated caracteristic times are 



r cool = Tdyn — mm (7jeans j Trot) > Tj eans = \ — — 7^—, T ro t = 7 T • (16) 

C V 327rGp v c {r) 

This formulation is usually termed the " Schmidt-Kennicutt prescription" . E and C are the 
gas internal energy and the effective cooling per unit volume (see equation [36]), while v c (r) 
is the galaxy rotational velocity at radius r. In the code the stars are maintained in the place 
where they form, and in each shell the associated sinks of momentum and internal energy 
per unit volume are given by the negative of 

• 4- ?7form^ ri+ VfoTmE /trr\ 

= , EJ = , (17) 

Tform Tf orm 
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where m is the specific momentum of the ISM (see Sect. 2.7). 

For a total mass AM* of newly formed stars in a given time-step and at a given place 
we assume a Salpeter IMF 

dN , ( M m [\ xl AM* ( M\~ x ~ x , , s 

so that the associated total number of Type II Supernovae is 

Nu = r d JL dM =(i--) (BA'^^Mi K 7 x 10 -3^, (19) 



Mu=8MQ dM V xJ\M a J M inf M M 

where the numerical value holds for x = 1.35. As for SNIa, we assume that each SNII event 
releases -Esn — 10 51 erg of kinetic energy, and the resulting mean efficiency is 



N u E SN r] SN /, 1\ ( MnfV M Q E SN rj SN 



'lais ^ i / -'-1111 i "dim '/din on., 1 n— 6 /on\ 

e n = — TTTT^- = 1 ~ Z HnH 77 — n^r - 3 - 9 x 10 Vj; (20) 



AM,c 2 V' xj \ M U J M inf M c 2 

in this paper we assume 7/sn = 0.85. The characteristic time for SNII explosion is fixed to 
Tn — 2 x 10 7 yr, and from equations (|15p and f[2"Uj) their luminosity (per unit volume) at each 
radius from the galaxy center is 

jfc(t) = — f pt{t')e- {t ' t,)/Tn dt'. (21) 



r n Jo 

We assume that each explosion leaves a neutron stars of 1.4M ; another possibility 
would be to assume more massive BH remnants (e.g., see Renzini & Ciotti 1993). As a 
consequence, the total mass ejected by the SNII explosions per unit mass is 

^(^r^.4^,0.2, (22) 
AM* \M n J AM* ' V ; 

and the mass return rate per unit volume of the young evolving stellar population is given 
by 

pn(t)~— f ptWe-W^M. (23) 
ni Jo 

Finally, in the code we also compute the fiducial optical and UV luminosity per unit 
volume of the new stars as 

E opt (t) = ^! f pt^e-^'y^dt', (24) 

Topt Jo 
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and 

Euv(t) = — f ptMe-W™**, (25) 

T"UV Jo 

respectively, where e op t = 1-24 x 10~ 3 , euv = 8.65 x 10~ 5 , T opt = 1.54 x 10 8 yr, and 
ruv = 2.57 x 10 6 yr are the efficiency and characteristic time of optical and UV emis- 
sion, respectively. The time-delay equations (1211) and (|23|) -( |25|) are integrated according to 
the scheme described in Appendix. 



2.4. The circumnuclear disk and the SMBH accretion luminosity 

At the onset of the cooling catastrophe a large amount of gas suddenly flows onto the 
central regions of the galaxy, and this induces star formation and accretion on the central 
SMBH, producing a burst of energy from the galaxy center. However observations of our own 
galactic center and high resolution studies of other nearby systems indicate that, in addition 
to the central starburst with radius ~ 100 — 300 pc, accretion onto the central SMBH is 
mediated by a small central disc within which additional significant star formation occurs 
(e.g., see Goodman & Tan 2004; Tan & Blackman 2005), and the remaining fraction of gas 
either is blown out in a broadline wind, or it is accreted onto the central SMBH. In our code 
the disk is not simulated with hydrodynamical equations, but its modelization is needed 
to obtain important quantities required by the code. The disk, which is the repository of 
the gas inflowing at a rate Mi from the first active grid point R\ of the hydrodynamical 
grid_|, and which feeds the central SMBH at a rate Mbh, contains at any time the mass gas 
m g and a total stellar mass m, = + m.*h> which is divided among low and high mass 
stars. Finally, the disk also contains a mass m rem of remnants from the earlier generations 
of evolved stars. 

In the adopted scheme the accretion rate on the central SMBH is given by 

Mbh = (26) 
1 + Vd 



where 



™fid = , ™Edd = Vd = 7T- > ( 27 ) 

r lag ec 2 2m Ed d 

are the fiducial depletion rate of gas from the circumnuclear disk and the Eddington mass 
accretion rate, respectively, while the characteristic disk star-formation time Ti ag is defined 



1 Mi is taken positive in case of accretion and zero in case of outflow at R\. 
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as 

2tt / Ft 3 , , s 

T '- s «imfc' (28) 

where a ~ 10~ 2 — 10 _1 (or higher, e.g., see King, Pringle & Livio 2007) is the disk viscosity 
coefficient. More rigorous formulations of ri ag are possible, but at the current level of mod- 
elization equation (28) is accurate enough. Thus, equations (26)-(27) guarantee that when 
rid << 1 the gas is accreted onto the central SMBH at the rate mgd, while Mbh = 2mEdd 
for r] d >> 1 (the factor of 1/2 in the definition of r] d allowing for possible super-Eddington 
accretion); note however that outside the first grid point R\ the flow accretion rate is limited 
in a self-consistent way by radiation pressure (see Sects. 2.6 and 2.7). 

We also assume that a fraction of the disk gas mass is converted into stars at a rate 
?7*^fid (where 77* ~ 10m g / Mbh> to approximate a Toomre criterion for star formation in the 
disk), and that another fraction of m g is lost as a wind at an instantaneous rate given by 
t] w Mbk (with r] w = 2, so that the broadline wind carries away twice as much mass as falls to 
the SMBH), so that the equation for the gas mass in the disk is 

dim 

— ± = Mi - (1 + r/ w )M BH - V*m& d . (29) 

The stars formed in the disk are described separately as a function of their mass, i.e., 
high-mass stars (M > Mu = 8M ) produce a total disk mass m*h, and low-mass stars 
(Minf < M < Mn) contribute to a disk mass m*i according to the equations 

—7— = (1 - hm* m ftd ; —7— = fhV^M , (30) 

at r*i at r*h 

where for the caracteristic evolutionary times we adopt r*i = r opt and r*h = rn given in 

Sect. 2.3, while we assume = 0.5, corresponding to a top-heavy Salpeter-like initial mass 

function of slope x ~ 1.16 and minimum mass M in f = 0.1M Q (see equation [IB]jE The 

associated optical (L^ opt ) and UV (Ld,uv) luminosities of the stellar disk are calculated 

following the scheme of Sect. 2.3. Finally stellar remnants mass in the disk evolves as 

dm iem _ m»i m*h 

1, Jrem,l ~r Jrem,h j \" / 

at r*i r*h 

where / rem ,i = 0.2, / rcmj h = 0.09. Thus, the equation for the mass of the disk wind is 

= ?7wM B H + (1 - /rem,l) h (1 ~ /rem,h) : (32) 

at r*i r^h 



2 The choice of a top-heavy IMF for the disk stars, at variance with the standard Salpeter law adopted for 
the new stars formed over the galaxy body, is motivated, for example, by the dynamical and X-ray evidences 
reported in Nayakshin & Sunyaev 2005 and Nayakshin et al. 2006). 
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the first term is a mass loss driven by the central SMBH and the second and third are from 
high mass and low mass stars in the central disk. All the equations in this Section are 
integrated with a first order finite difference scheme. 

From equation (26) we calculate the electromagnetic bolometric accretion luminosity as 

L BU {t) = eM BH (t)c 2 , (33) 

where the radiative efficiency e usually spans the range 0.001<e<0.15. We adopt e = 0.1 in 
the present paper, consistent with observations (e.g., see Soltan 1982, Yu & Tremaine 2002, 
Haiman et al. 2004), but a generalization to include an ADAF-like efficiency could also be 
easily implemented in the code (Narayan & Yi 1994, see also CO01). 

We also compute a fiducial mechanical luminosity for the disk wind Ldw, defined as 

£dw = eBHw^BHC 2 + e n c 2 (l - / rem ,h) — , (34) 

where e BHw = 5 x 10~ 4 = 0.005e is the mechanical efficiency of the SMBH (Elvis 2006). This 
mechanical energy output is a factor of 10 lower than adopted by Hernquist and collaborators 
when described in the same units, and adding this mechanical output to our simulation would 
be only a slight correction, since eeHw/e ~ 0.5%. The associated disk wind velocity is then 
given by 

2Ldw „ P^l c „ ? x 1Q 3 km B -i (35) 
m w V Vw 

for our parameters, in agreement with observations of broad line regions. 

In summary, the mass falling to the center is mediated by a gaseous, starforming a-disk 
which surrounds the SMBH. A small fraction of the mass is turned into (primarily) high 
mass stars in the disk, roughly two thirds is expelled in the broad line wind and roughly one 
third is accreted onto the central SMBH. 



2.5. Radiative heating and cooling 

With an improvement over CO01, the radiative heating and cooling produced by the 
accretion luminosity are numerically computed by using the Sazonov et al. (2005) formu- 
lae, which describe the net heating/cooling rate per unit volume E of a cosmic plasma in 
photoionization equilibrium with a radiation field characterized by the average quasar Spec- 
tral Energy Distribution by Sazonov et al. (2004), whose associated spectral temperature is 
T x ~ 2 keV. In particular, Compton heating and cooling, bremsstrahlung losses, line and 
recombination continuum heating and cooling, are taken into account. 
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A good approximation to the net gas energy change rate E, valid for T>10 4 K (all 
quantities are expressed in cgs system) is given by 



E = n 2 (S 1 + S 2 + S 3 ) = H-C, 



(36) 



where n is the Hydrogen density (in number), and positive and negative terms are grouped 
together in the heating (H) and cooling (C) functions. The bremsstrahlung losses are given 
by 

S 1 = -3.8 x 10~ 27 Vt, (37) 



the Compton heating and cooling is given by 

52 = 4.1 x 10~ 35 (T X — T) £, 



(38) 



where Tx is the Compton temperature, and finally the sum of photoionization heating, line 
and recombination continuum cooling is 



So = 10 
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(43) 



Equations (l38j) -(139l depend on the ionization parameter 



-^BH,photo( r ) 

n(r)r 2 



(44) 



where L^ u photo (r) is the effective accretion luminosity at r, which is evaluated by numerically 
solving in each shell the balance equation 



c ^BH,photo( r ) 

dr 



-Anr 2 H, 



(45) 
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with central boundary condition -^BH, P hoto( r = 0) = Lbh(£) given by equation fl33|) . The 
photoionization+Compton opacity associated with radiation absorption is then obtained 



Finally, the bolometric ISM luminosity is obtained from equation (I36I) as 

rr 

L r (r) = 4vr / CV 2 dr. (47) 



o 



The essential physics of this Section is well known. When the parameter £ is large 
thermodynamics guarantees that the gas temperature approaches the photoionization tem- 
perature ~ 10 7 ' 3 K, but for lower values of £ the temperature approaches ~ 10 4 K, near the 
peak of the cooling curve. 



2.6. Radiation pressure 

An important ingredient in the modelization of the gas flow evolution is the radiation 
pressure due to the accretion luminosity and to the light produced by the new stars. In 
its evaluation the explicit dependence on time is omitted, since the light travel time in the 
regions of interest is small compared to the time-scale on which the radiative input changes. 
Radiation pressure due to electron scattering (where neither the photon numbers, nor their 
energy change) is computed as 



K es pL BU + L vv (r) + Lopt(r) + L r (r) 

IVPradJes ~. n ' l 4i 

c 4nr z 



where n es = 0.35 in c.g.s. units, and from equations (j24|) - fl25]) 

L uv (r) = 4tt / E vv r 2 dr, L opt (r) = 4vr / E opt r 2 dr. (49) 
Jo Jo 

Note that all the luminosities used in equation (j4"8|) are unabsorbed. 
The radiation pressure contribution due to dust opacity is given by 

/V7 s KwP L BK,Vv( r ) + ^Uv( r ) K opt p^|H,opt( r ) + -^opt( r ) K m pL m (r) 
(Vp rad )dust = "j 2 a 2 A 2~' ( 50 ) 

c Airr z c 47rr z c Anr z 



where 



L IR (r) ee L^ )UV (r) + L^ >opt (r) + Lg»(r) + L*», (51) 
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is the infrared luminosity due to recycling of photons absorbed from the ISM, and we adopt 
as estimates for (cgs) opacity in three bands 

K o P t = J — jy-j_Q4' K uv = 4/t pt, Kir = ^g, (52) 

where the temperature dependent denominator is designed to mimic the destruction of dust 
at high temperatures (T. Thompson & B. Draine, private communication); the dust opacity 
we are using is likely a lower bound to the true value, while a more accurate treatment 
can be found in Thompson, Quataert & Murray (2005), and will be implemented in future 
explorations. At variance with electron scattering the effective luminosities appearing in 
equations ( |50l - (l5B take into account absorption, and are obtained by numerically solving 
the two lowest spherically symmetric moment equations of radiative transfer in the Eddington 
approximation (e.g., Chandrasekhar 1960): 



J T eff dL eS 

= 47rr 2 (E uv - KwpJfy), = 4vrr 2 (£ opt - n optP J^ t ) 



dJ\jy _ _3^uvA^uv ^opt _ 3K ptpLp pt 
dr 47rr 2 ' dr Airr 2 ' 

The central boundary conditions for stellar luminosities are L^(0) = Ld,uv; -^opt(O) = £d,opt> 
jeff (o) = L^uv/W^Rl and Jg(0) = L A ^ t /\QTX 2 R\ (see Sect. 2.4). The effective accretion 
luminosities -^bh,uv an d -^BH, op t are computed with two equations similar to (1531) . where 
the distributed source term is missing, J = L^/inr 2 , and in the UV and optical bands 
£bh,uv(°) = °- 2L B H (t) and Lf UjOpt (0) = 0.1L BH (t), respectively. 

The last contribution to radiation pressure comes from photoionization opacity, 



reft 

P^photo -^BH 



r 



I VPradJphoto — : n , (.00) 

c 4irr z 

where the photoionization opacity and the absorbed accretion luminosity are calculated as 
described in Sect. 2.5. It is well known that the radiation pressure on electrons (equation [45] ) 
can significantly affect the gas dynamics when the AGN luminosity approaches the Eddington 
limit. Consistently with the findings of Thompson, Quataert, & Waxman (2006) we also find 
that the radiation pressure on the dust can have a major effect during starburst phases in 
retarding the infall of cool gas, thus boosting the mass transformed into stars and reducing 
the gas available for accretion onto the central SMBH. 



2.7. Hydrodynamical equations 



As in CO01, the evolution of the galactic gas flows is obtained integrating the time- 
dependent Eulerian equations of hydrodynamics, where now we have several additional source 
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and sink terms 



dm 
~dt 



dp 

— + V ■ (pv) = ap* + pn- p+, (56) 
V ■ (mv) = -(7 - 1)VE - Vp rad + gp- m+ (57) 



<9£ 

— + V-(£to) = -(7-l)SV-u+ff-C+ (58) 



(ap* +Pii)(y 2 + 3<?1) , p , p ^ 
h -C/i + -C/ii — 



p, m, and are the gas mass, momentum and internal energy per unit volume, respectively, 
and v is the gas velocity. The ratio of the specific heats is 7 = 5/3, and g(r) is the grav- 
itational field of the galaxy (stars and dark matter), plus the contribution of the central 
SMBH. The gravitational field is updated at each time step by considering the SMBH mass 
growth; for simplicity, we do not take into account neither the ISM contribution, nor the 
mass redistribution due to the stellar mass losses and star formation. The total radiative 
pressure gradient is Vp rad = (Vp rad ) e s + (Vp rad ) dust + (Vp rad ) ph oto (Sect. 2.6), while the 
radiative heating and cooling term H — C is described in Sect. 2.5. 

The energy source term is obtained under the assumption that the streaming velocity 
of the source distribution is zero, neglecting the small contributions of the internal energy 
of the injected gas, and of the kinetic energy of stellar wind when compared to the local 
stellar velocity dispersion contribution (for the derivation and detailed discussion of the 
hydrodynamical equations with moving isotropic or anisotropic source terms, see Recchi, 
D'Ercole & Ciotti 2000). Note that the term proportional to the stellar velocity dispersion 
becomes dominant near the SMBH, as described in equation (J7|). The source terms ap* and 
Ej of the initial, passively evolving stellar population are given in equations (jT31 -( fT4l) . while 
the source terms due to Type II Supernovae, p\\ and En, are described in Sect. 2.3. 

From the numerical point of view, the code is the same as in C097 and CO01, however 
the present simulations are much more difficult to run, because the stellar density distribu- 
tion, being more concentrated than the King model, now injects more gas into the central 
regions, thus increasing the gas density and decreasing its cooling time. In addition, the 
gas at the galaxy center is gravitationally more bound (due to higher mass concentration of 
the new models), and correspondingly more difficult to expell. Finally, the numerical grid 
spacing has been reduced in order to resolve the central regions of the galaxy. In particular, 
we place the first active grid point R\ within the Compton radius 

2GM BH pm p M BH 10 7 K 

Rx = 3k B T x " 3 - 6/i 10^^T PC ' (59) 
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so that at R\ we can impose the physical condition of a vanishing thermodynamical pressure 
gradient, leading to gas free-fall on the circumnuclear disk when the radiation pressure is 
negligible; in this paper we adopt Tx = 2.5 x 10 7 K. The appropriate values for radiation 
pressure at R\ are obtained from the disk treatment described in Sects. 2.4 and 2.6. Note 
that in CO01 we were not able to perform a full simulation with such high resolution, as a 
consequence of the higher Tx adopted. 

The simulations are realized with a spatially second-order Eulerian scheme which adopts 
two staggered grids (for scalar and vector quantities, see CDPR and CO01 for details), each 
of them consisting of 120 logarithmically spaced grid points, covering the range 2.5 pc - 
200 kpc. The equations are integrated with a time-splitting scheme, while the heating and 
cooling terms in the energy equation are integrated by using a predictor-corrector scheme, 
so that the integration is second order in time. At each simulation time, the time-step is 
determined as a fraction of the minimum among the Courant condition over the grid, and 
of the others characteristic times associated with the described physical processes: during 
the accretion phases (and subsequents bursts of radiation), it is not infrequent to have time- 
steps of the order of 1 yr or less. However, it is important to note the accretion events are 
characterized by the intrinsic time-scale related to equation (1511 by 

where cx is the isothermal sound velocity associated with the Compton temperature. 

3. Model evolution 

We now show the main properties of a representative model characterized by an initial 
stellar mass M* = 4.6 x 10 11 M Q , a FP effective radius R e = 6.9 kpc and aperture velocity 
dispersion cr a = 260 km s _1 (leading to an expanded initial condition of R c = 10.35 kpc and 
o~ a = 235 km s _1 ), total dark-to- visible mass ratio TZ = 5 and dark-to- visible scale- length 
ratio j3 = 5.22 (corresponding to an identical amount of stellar and dark matter within 
the half-light radius). The initial SMBH mass follows the present day Magorrian relation, 
i.e., Mbh — 10~ 3 M*. We remark again that this model galaxy is not fully appropriate 
as an initial condition for a cosmological simulation, because its parameters are fixed to 
reproduce an early-type galaxy similar to those observed in the local universe, and also 
because we set outflow boundary conditions at the galaxy outskirsts (~ 200 kpc): from this 
point of view, the simulations represent an isolated elliptical galaxy (note for example that 
in the present context we are not considering the effects of possible merging on the galaxy 
evolution). We adopted this procedure to adhere to the standard approach followed in 
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" cooling-flow" simulations, while in future explorations we will address in a more consistent 
way the problem of the galaxy structural and dynamical modifications due to star formation 
and mass redistribution over an Hubble time, and the compatibility of the obtained galaxies 
with the present-day scaling laws of elliptical galaxies. We remark that the model presented 
is just one out several tens of runs that have been made, characterized by different choices 
of the parameters (often outside the currently accepted ranges), for example with high/null 
DM, enhanced/suppressed star formation, vanishing efficiencies, and so on. The obtained 
model evolution can be very different, some of them leading to galaxies in a permanent 
wind state, or galaxies with extremely massive final SMBHs. A discussion of these issues is 
deferred to successive works. 

The initial conditions are represented by a very low density gas at the local virial 
temperature. The establishment of such high-temperature gas phase at early cosmological 
times is believed to be due to a "phase-transition" when, as a consequence of star formation, 
the gas-to-stars mass ratio was of the order of 10% and the combined effect of SNIa explosions 
and AGN feedback became effective in heating the gas and driving galactic winds. Several 
theoretical arguments and much empirical evidence, such as galaxy evolutionary models and 
the metal content of the Intra Cluster Medium (ICM) support this scenario (e.g., Renzini 
et al. 1993; OC05; Di Matteo et al. 2005). For the reasons above, in the simulation here 
presented (as well as in all others simulations not shown), we assume that the age of the 
galaxy stellar component at the beginning of the simulation is 2 Gyr old, and the simulations 
span 12 Gyr, so that the cosmic time at the end of the simulations is 14 Gyr. 



3.1. Luminosities 

A first, important result of the new models is that overall the main properties of the 
CO01 models are confirmed. After a first evolutionary phase in which a galactic wind is 
sustained by the combined heating of SNIa and thermalization of stellar velocity dispersion, 
the central "cooling catastrophe" commences. In absence of the central SMBH a "mini- 
inflow" would be then established, with the flow stagnation radius (i.e., the radius at which 
the flow velocity is zero) of the order of a few hundred pc to a few kpc. These "decoupled" 
flows are a specific feature of cuspy galaxy models with moderate SNIa heating (Pellegrini 
& Ciotti 1998). However, after the central cooling catastrophe, the feedback caused by 
photoionization and Compton heating strongly affects the subsequent evolution, as can be 
seen in Fig. [T]where we show the luminosity evolution of the central AGN with time-sampling 
of 10 5 yrs. The bolometric luminosity (top panel) ranges between roughly 0.1 to 0.001 of the 
Eddington limit (the almost horizontal solid line) at peaks in the SMBH output but, since 
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Fig. 1. — Dotted lines are the bolometric accretion luminosity (top panel) and the optical 
SMBH luminosity corrected for absorption -^bh opt > i- e > as would be observed from infinity 
(bottom panel, see Sect. 2.6); we recall that at the center we fixed L^ H t (Ri) = 0.1L B h- 
The almost horizontal solid line is L Edd . 
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Fig. 2. — Panel a: time expansion of 500 Myrs of top panel in Fig. [H showing the first major 
burst. The dotted line is the bolometric accretion luminosity. The widths of the spikes is 
typically of the order of ~ 0.1 — 1 Myr. Total opacities (defined as r = J Q Rt npdr, where 
R t = 200 kpc) of dust on UV (dotted) and optical (dashed) luminosities (panel b), and 
electron scattering (solid), and photoionization opacity (long dashed line) (panel c). Dust 
opacity on the IR recycled radiation would be a line parallel to those in panel b, but ~ two 
order of magnitude lower (see equation [52]). Note that the vertical scale is the same in 
panels b and c. 
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obscuration is often significant, the optical accretion luminosity as seen from infinity can be 
much lower still (bottom panel). But the central quasar is not always obscured and we see, 
in the lower panel of Fig. [TJ that the optical luminosity reaches ~ 10 44 erg s -1 in numerous 
bursts. As already found in CO01, the major AGN outbursts are separated by increasing 
intervals of time (set by the cooling time), and present a characteristic temporal substructure, 
whose origin is due to the cooperating effect of direct and reflected shock waves. At t ~ 8 
Gyrs the SNIa heating, also sustained by a last strong AGN burst, becomes dominant, a 
global galactic wind takes place and the nuclear accretion switches to the optically thin 
regime. Note that a further reduction of the accretion luminosity during this phase would 
be otained if considering ADAF accretion instead of standard accretion. 

The temporal substructure of the first major burst is revealed in Fig. [2h, where we show 
a blow-up of 500 Myr of top panel Fig. [Tj, starting at 2.5 Gyr: the time extent of each of the 
sub-bursts (for example, when Lbh > 10 45 erg s _1 ) is of the order of ^1 Myr. In Fig. [2b, c we 
also show the time evolution of the different model opacities during the same burst, where 
the transition from optically thin to the optically thick and back to thin phase is apparent. 

In Fig. [3h we show the coronal X-ray luminosity Lx (emitted by gas at T > 5 x 10 6 
K), due to the hot galactic atmosphere integrated within 10i? e . In particular, Lx falls in 
the range commonly observed in massive early-type galaxies, with mean values lower than 
the expected luminosity for a standard cooling-flow model. In Fig. [3b we show instead the 
estimated IR luminosity L IR due to the reprocessing of the radiation emitted by the new stars 
and by the SMBH and absorbed by the ISM inside 10R e (equation [51]). The simulations 
show that the bulk of the reprocessed radiation comes from AGN obscuration, while the 
lower envelope is set by radiative reprocessing from the new stars. The very high luminosity 
peaks (L IR ~ iq45-46 er g g -i^ correS p 0nc l to one component of the SCUBA sources seen at 
z ~ 2 (e.g., see Pope et al. 2006). Finally, in Fig. [3fc we show the temporal evolution of the 
optical and UV luminosities of the starbursts corrected for absorption. Overall, Figs. [3bc 
show that a large fraction of the starburst luminosity output occurs during phases when 
shrouding by dust is significant (e.g., see Rodighiero et al. 2007), i.e. the model would be 
observed as an IR source with UV and optical in the range seen in brighter E+A sources. 

The spatial radius within which half of the IR and X-ray luminosities are emitted changes 
dramatically with time, and as a function of the total emitted luminosity. This is apparent 
from Fig. HJ where we show the time evolution of r^x and rh,m. Note the large variation 
of the size of the emission regions: in particular, the smallest values of the half-light radii 
correspond to the peaks of the associated luminosities. During strong radiation outbursts, 
the emitting regions are so small (~ 10 pc or less) that virtually the entire X-ray and IR 
luminosities would be seen as emitted by the central source; these phases also correspond to 
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Fig. 3. — The galaxy X-ray coronal luminosity Lx (top), recycled infrared Lir (middle), and 
the starburst UV and optical luminosities (bottom), corrected for absorption. 
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phases of high obscuration, with optical depth of the order of 100 (e.g., Imanishi et al. 2007). 
Instead, during the late-time hot accretion phase r^x reaches values commonly observed in 
elliptical galaxies, and -Rh,m is even larger. When the IR luminosity is as large as seen in 
the SCUBA observations (Pope et al. 2006), we predict that the characteristic sizes will be 
of the order of ~ 10 2,5 pc. 

An important quantity associated with the time evolution of the various luminosities is 
their duty cycle. As in CO01 for a given luminosity L(t), we define the associated duty-cycle^] 
over a period of time At 

h "> = (61) 

In Fig. [5]we show the resulting values for a time dependent temporal window At = t/2. 
In the top panel we show the (logarithmic) value for the effective optical and UV AGN 
luminosities (with negative peaks corresponding to values as low as 0.01 or less), while 
during the smooth accretion of the last Gyrs the values flatten to unity. In the bottom 
panel the duty-cycles correspond to the starburst optical and UV luminosities, and show 
a larger and less fluctuating values ^0.5, in agreement with observational results (e.g., see 
Cimatti et al. 2002). In the same panel we also show the duty-cycle of the X-ray ISM 
luminosity computed outside a sphere of radius 100 pc, so as to exclude the ISM luminosity 
fluctuations produced by direct AGN heating on the surrouding ISM. In other words, the 
X-ray duty cycle refers to the bulk of the galaxy body, and should give an indication of the 
expected fraction of significantly disturbed galaxies in coronal X-rays. Quite obviously, the 
derived values depend on the adopted sampling time interval At. For example, by fixing 
the sampling time to A = 100 Myrs, the values in the bottom panel are almost unchanged, 
while the AGN duty-cycles becomes (during peaks of activity) as small as ~ 10~ 3 . This is 
consistent with the temporal substructure of the major bursts (see Fig. [2]) At the opposite 
end, taking all the time interval spanned by the simulation, the AGN duty-cycle (both in 
UV and optical) is ~ 10~ 2 , the IR is ~ 0.2, while the starburst duty-cycle is ~ 0.8 and that 
of the global ISM X-ray ~ 0.4. We stress here that the duty-cycles are computed by the 
code calculating the luminosity values at each time-step, that is usually of the order of the 
year or even less. 

Duty cycles can be also defined in a different way. For example, in Table 1 we focus 
in particular on the estimated time fraction spent by the central SMBH shining at a given 



3 Defined in this way fduty would be the fraction of the time a system would be in a high-luminosity 
state, L/i, and 1 — fduty would be the fraction of the time it spent in a low-luminosity state Li, if it were to 
oscillate between these two states and Li/Lh « 1. 
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Fig. 4. — Time evolution of the (volume) half-light radius of the X-ray ISM luminosity 
(bottom panel) and IR luminosity (top panel) during the bursting phases. Small radii 
correspond to the high-luminosity peaks, and the predicted SCUBA-like sources should be 
of linear size ~ 10 2 5 pc. 
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Fig. 5. — Time evolution of duty cycles as computed from equation (61), with a time-window 
At = t/2. Top panel: duty-cycle of £bh,uv (solid) and L^ opt (dotted); the top axis shows 
the corresponding redshift (e.g., Spergel et al. 2006). We see that these systems would 
be observed from afar in the (rest-frame) optical or UV as quasars several percent of the 
time. Bottom panel: duty-cycle of the starburst L^f Y (solid), L^ t (dotted), of the ISM X- 
ray luminosity (computed in a volume excluding the inner 100 pc), and of the recycled IR 
luminosity L IR . 
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fraction of Eddington luminosity, also considering the obscured accretion phases. Overall, 
in the bursting phase (1<£<3), the duty-cycle of the SMBH in its "on" phase is of order 
precents and it is unobscured approximately one-third of the time. We found it interesting 
that these figures, obtained from hydrodynamical simulations, can be positively compared 
with observations (e.g., see Gilli, Comastri & Hasinger 2007; Martinez-Sansigre & Rawlings 
2007). 

3.2. Mass budgets 

In Fig. Owe show the time evolution of some of the relevant mass budgets of the model, 
both as time-integrated properties and instantaneous rates. At the end of the simulation 
the total ISM mass in the galaxy is ~ 5 x 1O 8 M , while the SMBH mass increased by 
2.5 x 1O 8 M , thus reaching a final mass of ~ 7 x 1O 8 M : a model with a smaller initial 
SMBH mass would accrete less, thus maintaining the Magorrian relation even better. The 
SMBH mass accretion rate strongly oscillates as a consequence of radiative feedback, with 
peaks of the order of 10 or (more) M /yr, while during the final, hot-accretion phase the 
almost stationary accretion is ^10 _4 M Q /yr: this value is close to the estimates obtained 
for the nuclei of nearby galaxies (Pellegrini 2005). Note that in the last 6 Gyrs the SMBH 
virtually stops its growth, while the ISM mass first increased due to the high mass return 
rate of the evolving stellar population, and then decreases due to the global galactic wind 
induced by SNIa. During the entire model evolution, >lO ia5 M of recycled gas has been 
added to the ISM from stellar mass losses. Approximately 2.1 x 1O 1O M has been expelled 
as a galactic wind, while ~ 1.4 x 1O 1O M has been transformed into new stars, so that only 
0.7% of the recycled gas has been accreted onto the central SMBH, and the central paradox 
of the mass budget is automatically resolved. 

In the present simulation, approximately twice of the total mass accreted onto the central 
SMBH is expelled as a disk wind, while the final mass of the disk, in stellar remnants, sums 
up to m Icm ~ 2.8 x 1O 5 M . It is important to stress that an identical model without SMBH 
feedback (i.e., e = in equation [33]), but with the same star formation treatment of the 
model described in this paper, produced a SMBH of final mass >1O 1O M , while the total 
mass in new stars was reduced to ~ 3 x 1O 9 M . In addition, this "ad hoc" model does 
not present fluctuations in the starburst and ISM X-ray luminosities, thus showing the vital 
importance of SMBH feedback in the overall results. Tests of numerical convergence were 
performed to determine the extent to which the results quoted in this paper would be altered 
as one increased the spatial and temporal resolution. To this end we reduced the grid spacing 
by a factor of 1.5 and then a factor of 2.0. Almost all results changed at the level of a few 
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Fig. 6. — Mass budget evolution. Top left panel: total hot gas mass in the galaxy (within 
10R e , M ISM , dotted line), and accreted mass on the central SMBH (AM BH ). Top right: 
mass lost as a galactic wind at 10R e (AM^, dashed line), and total mass of new stars 
(AM*) formed according to equation (fT5|) . Bottom left panel: mass return rate from the 
evolving stellar population (as given by volume integral of equation [15] . dotted line), and 
mass accretion rate on the central SMBH (Mbhj equation [26]). Bottom right: galactic 
wind mass loss rate at 10R e (M^, dashed line), and instantaneous, volume integrated, star 
formation rate. Note that for t > 8 Gyrs the mass lost as a galactic wind is almost coincident 
with the mass input from evolving stars. 
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percent or less (with the numerical uncertainty thus far below the level of the uncertainty in 
the physical modeling). The largest change found in the highest resolution run was in the 
growth of the central SMBH mass with this growth being reduced by ~ 14% . This is in 
the direction to make the final results conform more closely to the Magorrian relation with 
regard to AM BH /AM t (see Fig. (Bj). 

Two important quantities associated with the mass budget of the model are the me- 
chanical (non-relativistic) wind luminosity of the disk as given by equation (1341) . and the 
mechanical luminosity of the galaxy, i.e. the kinetic energy carried away by the global galac- 
tic wind (Fig. 0, bottom panels, while the corresponding mass rates are shown in the top 
panels). Time integration of these two mechanical luminosities over the entire model evo- 
lution revealed that the disk wind would deposit ~ 2.2 x 10 59 erg in the galaxy ISM, while 
the galactic wind would inject in the ICM ~ 1.3 x 10 58 erg. The star formation rate during 
the periods of feedback dominated accretion oscillates from 0.1 up to several hundreds (with 
peaks near 10 3 ) M® yr -1 , while it drops monotonically from 10 _1 to ^10~ 3 M© yr _1 in the 
last 6 Gyrs of quiescent accretion (see Fig. 6). As already mentioned above, these violent star 
formation episodes (with SMBH accretion to star formation mass ratios ~ 10~ 2 or less) are 
induced by accretion feedback^, and are spatially limited to the central 10 — 100 pc; thus, 
the bulk of gas flowing to the center is consumed in the starburst. These findings are nicely 
supported by recent observations (e.g, see Sect. 5 in Lauer et al. 2005, see also Davies et 
al. 2007). Note that the "age" effect of the new stars on the global stellar population of the 
galaxy is small, as the new mass is only 3% of the original stellar mass, and it is virtually 
accumulated during the first Gyrs (see Fig. [6]), so that the mass- weighted age of the final 
model is of the order of 12 Gyrs. The half-mass radius of the final stellar distribution (with- 
out considering adiabatic contraction, nor the reduction of the stellar mass distribution due 
to galactic winds, see Section 4) contracts by ~ 16%, just due to the addition of the new 
stars in the central regions of the galaxy. This is made apparent in Fig. [SJ where we show 
the final spatial density profile of the system, together with its projection and the best-fit 
obtained with the Sersic (1968) law 

S( J R) = S e- b ( m )^) 1/m , (62) 

where b(m) = 2m — 1/3 + 4/(405m) (Ciotti & Bertin 1999). As expected, the profiles show 
an increase of the best-fit Sersic parameter m, due to the mass accumulation in the central 
regions. Remarkably, the final value of m is within the range of values commonly observed 
in ellipticals: however, in the final model we note the presence of a central (~ 30 pc) nucleus 



4 However, bursting star formation is not necessarily associated with AGN feedback (e.g., see Kriigel & 
Tutukov 1993). 
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Fig. 7. — Top panels: temporal evolution of the disk wind mass rate (equation [32], left) 
and the galactic wind mass rate computed at lOR e , M gw = 4n(10R e ) 2 p(10R e )v(10R e ) (right). 
Bottom panel: the corresponding mechanical luminosities, i.e., the kinetic energy that would 
be injected from the disk in the galaxy interstellar medium (L^ w ), and by the galaxy in the 
ICM, = M gw (10 J R e )t; 2 (10/? c )/2. 
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Fig. 8. — Dotted lines are the projected surface density of the model shortly after the 
beginning of the simulation (t = 2.5 Gyr, z ~ 2.6, bottom panel), and at t — 13.5 Gyr 
[z — 0, top panel), normalized to the surface density at the effective radius. The relation 
between age and redshift holds for standard Cosmology (e.g., Spergel et al. 2006). Solid 
lines are the best-fit Sersic law. The effective radius contracted from ~ 9.2 kpc to ~ 8.4 kpc, 
while the surface density S e increased from ~ 3 x 10 22 to ~ 3.6 x 10 22 protons per cm~ 2 . 
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originated by star formation which stays above the best fit profile, similar to the light spikes 
characterizing "nucleated" galaxies (e.g., see Graham & Driver 2005, Lauer et al. 2005). 

3.3. Hydrodynamics 

In Fig. Owe show the temperature and density in the central regions of the model: note 
how the SMBH bursts heat the central gas, causing the density to drop, and launching gas at 
positive velocities of the order of thousands km s _1 (this can better appreciated in Fig. [TOj 
where we show a time blow-up of the first two SMBH feedback events). The Compton 
temperature is the horizontal dashed line, and during the bursts the local gas is heated 
above this limit. 

As was already found in CO01, the galaxy cooling catastrophe starts with the formation 
of a cold shell placed around the galaxy core radius: however, in the present models (as in 
those explored by Pellegrini & Ciotti 1998), the cooling catastrophe happens at significantly 
earlier times than in CO01 models, both due to the higher central stellar density and to 
the different time dependence and amount of SNIa explosions. The three main evolutionary 
phases of the model are summarized in Figs. fTTlTbil In particular, in Fig. [11] one can observe 
(with time separation of 1 Myr), the evolution of the first cold shell falling to the galaxy 
center, while in Fig. [T2] (dashed line) the expanding material due to the first burst (with 
velocities of the order of several hundred km s _1 ) is clearly visible. A particularly important 
feature can be noticed in Fig. [131 where a new cold and dense shell of gas is formed as a 
consequence of the shock wave produced by the burst. This shell moves (slowly) outwards, 
and then starts to fall back at the center; in the shell, star formation rate reaches values of 
~ 10 _8 M Q yr _1 pc -3 . We stress that this shell is of a different origin compared to the first, 
and it is possibly Rayleigh- Taylor unstable. This cycle of shell formation, central burst, 
and expanding phase, repeats during all the bursting evolution, along the lines described 
in detail in CO01. Finally, when the specific SNIa heating becomes dominant over the 
decline of fresh mass input from evolving stars, the galaxy hosts a wind, and the accretion 
becomes stationary without oscillations, and the central SMBH is radiating at ~ 10~ 5 LEdd 
(Hopkins, Narayan & Hernquist 2006). The hydrodynamical quantities (separated now by a 
time interval of 0.5 Gyrs) are shown in Fig. [14j 
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Fig. 9. — Top panel: gas velocity at 5 pc from the SMBH. Note how the SMBH growth affects 
the lower envelope of velocity values. Bottom panel: Gas number density (dotted line, scale 
on the right axis) and temperature at 5 pc from the SMBH (solid line). Low-temperature, 
high-density phases end when accretion luminosity Lbh increases sharply heating the ambi- 
ent gas to a high-temperature, low-density state. The horizontal dashed line is the model 
Compton temperature T x = 2.5 x 10 7 K. 
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Fig. 10. — Time expansion of the first two feedback events of the initial major burst in Fig. [9] 
(whose total time extent is ~ 400 Myr) at time resolution 10 5 yr; the horizontal dashed line 
is the model Compton temperature. The complementary behavior of p and T is apparent. 
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Fig. 11. — Hydrodynamical quantities of the model immediately before the onset of the 
bursting phase (2.554 Gyr), separated by 1 Myr (in order, solid, dotted, dahsed). The cold 
shell is falling towards the galaxy center. 
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Fig. 12. — The cold shell reached the center, and a shock wave is moving outward (dashed 
lines, note the positive velocities). Time interval is 1 Myr (in order, solid, dotted and dashed 
lines), and the first two snapshots are the two last snapshots in Fig. [TTJ 
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Fig. 13. — The expanding flow produced a new cold shell, which starts to fall to the galactic 
center. Time interval is 1 Myr (in order, solid, dotted and dashed lines), while the time is 
now 2.561 Gyr. 
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Fig. 14. — Hydrodynamical quantities of the model at the end of the simulation, separated 
by 0.5 Gyrs (in order, solid, dotted and dashed lines). The galaxy is in the low-luminosity, 
high-temperature and low-density stationary accretion phase in the inner parts and, in its 
outer parts a ~100-200 km s _1 wind is carrying out most of the late recycled gas. 
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4. Conclusions and discussion 

In this paper we have addressed the effects of radiative feedback on the gas flows in 
elliptical galaxies. The investigation is in the line of previous exploratory papers (C097, 
CO01, OC05), but now the input physics and the galaxy models have been substantially 
improved. We briefly recall here the main points on which our framework is based. First 
of all, it is obvious that the recycled gas from dying stars is an important source of fuel for 
the central SMBH, even in absence of external phenomena such as galaxy merging, that are 
often advocated as the way to induce QSO activity. It is also obvious that the recycled gas, 
arising from stars in the inner several kpc of the galaxy (assumed a giant elliptical), will 
necessarily drive a classical radiative instability and a collapse towards the center of metal 
rich gas. As a consequence, a star-burst must occur and also the central SMBH must be 
fed. The details of how much is accreted on the central SMBH vs. consumed in stars vs. 
ejected from the center by energy input from the starburst and AGN are uncertain. But 
order of magnitude estimates would have the bulk going into stars or blown out as a galactic 
wind, with a small amount going into the central SMBH. In addition, since at the end of 
a major outburst an hot bubble remains at the center, both processes shut themselves off, 
and it will take a cooling time for the cycle to repeat. In other words, relaxation oscillations 
are to be expected, but their detailed character is uncertain. Finally, order of magnitude 
estimates would indicate that during the bursting phase the center would be optically thick 
to dust, so one would observe a largely obscured starburst and largely obscured AGN with 
most radiation in the far IR; as gas gets consumed, the central sources become visible. Much 
of the AGN output occurs during obscured phases: then there is a brief interval when one 
sees a " normal" quasar, and finally one would see a low X-ray luminosity and E+ A spectrum 
galaxy, with A dominating in the central several hundred pc for 10 7-8 yrs. 

The present paper attempts to illustrate the expectations described above. Overall, we 
have confirmed that radiative feedback from a central SMBH has dramatic effects on galaxy 
evolution and on the mass growth of the central SMBH itself, and we find that much of the 
recycled gas falling towards the galaxy center during the accretion events is consumed in 
central starbursts with a small fraction (of the order of 1% or less) accreted onto the central 
SMBH. In particular, in the presented model approximately half of the recycled mass from 
evolving stars is expelled in the ICM, while the remaining fraction is consumed in central 
starbursts. Thus, the central starburst is an important component in the physical modelling, 
since it regulates the amount of gas available to be accreted onto the central SMBH: without 
allowing for the (AGN feedback induced) central star formation the SMBH would grow to be 
far more massive than seen in real galaxies. While the details predicted by our simulations 
are uncertain, we do not see how this sequence can be avoided in its qualitative features. The 
input is standard physics plus the known radiative output from the SMBH and stars. Other 
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processes that we do not include in the present code may also be important. For example, 
the mechanical energy from the AGN and cosmic rays from SNR in the starburst is surely 
important (e.g., see Di Matteo et al. 2005), but we deliberately excluded them from the 
present investigation to better assess the importance of radiative feedback (both as heating 
and radiative pressure). 

The main results of our simulations, also considering all the simplifications in the treat- 
ment of physics, and of the geometry of the code, nicely support the scenario depicted above. 
In particular, we showed how complicated the evolution of an isolated galaxy, subject to inter- 
nal evolution only, can be (e.g., see Pierce et al. 2007). From the observational point of view, 
we find that evidence for starbursts should be common when looking at elliptical galaxies, 
with the fraction showing an E+A spectrum increasing with redshift (with preliminary evi- 
dence for high metallicity in the new stellar spectra). Also, signs of star- formation in high-z 
objects such as those detected in UV (GALEX) and IR (Spitzer), or accompanying AGN 
(X-ray), should be common (e.g., see Yan et al. 2006; Nesvadba et al. 2006, Simoes Lopes et 
al. 2007), with IR luminosity peaks of ~ 10 46 erg s -1 . Due to the observational relevance of 
these predictions, a more accurate modelization of the expected IR model properties would 
be a natural extension of the present work (e.g., Chakrabarti et al. 2007). 

We note that there is increasing evidence in the local universe of hot gas disturbances 
on various galactic scales, likely the resultant from recent nuclear activity (e.g., Forman et 
al. 2006). For example, Chandra revealed two symmetric arm-like features across the center 
of NGC4636 (Jones et al. 2002; O'Sullivan, Vrtilek, & Kempner 2005), accompanied by a 
temperature increase with respect to the surrounding hot ISM; they were related to shock 
heating of the ISM, caused by a recent nuclear outburst. Other evidences include a hot 
filament in the nuclear regions of NGC821 and NGC3377 (Fabbiano et al. 2004, Soria et 
al. 2006); a "bar" feature, presumably due to a shock, at the center of NGC4649 (Randall, 
Sarazin & Irwin 2004); a nuclear outflow in NGC4438 (Machacek, Jones, & Forman 2004); 
an unusual temperature profile finally is present in NGC 3411 (Vrtilek et al. 2006). 

In this paper we presented just one out many models, to illustrate in detail the global 
behavior of a typical solution. Of course, the present work suffers of some weak points, as: 
1) we neglected the modifications of the galaxy gravity field and velocity dispersion profile 
due to the stellar mass losses, galactic wind, and star formation (deepening the central 
potential). 2) The new stars are placed in the galaxy where they form. 3) The additional 
effect of mechanical feedback was not taken into account (e.g., see Begelman & Ruszkowski 
2005). 4) Finally, the simulations are spherically symmetric, so that Rayleigh- Taylor unstable 
configurations of the ISM, and the formation of an accretion disk, cannot be followed. 

Several lines of investigation will be studied in future papers, in order to better test the 



-40- 



results and to address specific points only mentioned here. For example, it is natural to study 
in detail the properties expected for the starburst population (such as spatial distribution, 
spectral properties, etc.), and the X-ray properties of the perturbed ISM, as a function of the 
combined effect of SNIa and central feedback. Other obvious issues to be addressed are the 
effects of environment, as for a cD galaxy in a cluster; this study will also probe the impact of 
photoionization plus Compton heating on the ICM (extending the preliminary investigation 
of Ciotti, Ostriker, & Pellegrini 2004, see also Pope et al. 2007). Also, the modifications of 
the galaxy inner structure and dynamics due to star formation are important, in particular 
for the expected evolution of the Magorrian and Mbr—u relations. A most relevant extension 
of the present work would also be the study of radiative feedback by using two-dimensional 
codes, in order to follow the evolution of unstable features here found, such as the cold-shell 
phase, and to properly describe axisymmetric accretion in the optically thick regime (e.g., 
Krolik 2007, Proga 2007). We finally mention another observational riddle that could be 
addressed with the present models, i.e., that of the apparent "underluminosity" of SMBHs 
in the local universe (e.g., see Fabian & Canizares 1988; Pellegrini 2005). The present 
models could also be used in a different way. In fact, when addressing galaxy evolution 
over cosmo logical times, two main effects have not been considered in the present treatment, 
namely galaxy formation itself and the successive cosmological infall and merging events 
(e.g., Hopkins & Hernquist 2006, Micic et al. 2007). These aspects of evolution are usually 
considered in cosmological simulations (e.g., see Springel, Di Matteo, & Hernquist 2005; 
Naab et al. 2007, and references therein), and so it would be interesting to add at the center 
of galaxies in those simulations the treatment of feedback physics here described. 

We thank Annibale D'Ercole, Bruce Draine, Martin Elvis, Roberto Gilli, Jeremy Good- 
man, Avi Loeb, Eve Ostriker, Silvia Pellegrini and Todd Thompson for useful discussions, 
and the anonymous Referee for important remarks. L.C. was partially supported by the 
grant CoFin2004 by Italian MIUR. 
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A. Numerical evaluation of lag integrals 

Time-lag integrals often appear in the code. It is possible to compute them numerically 
without much use of computer memor}@. In fact, let f(t) a given function of time, and let 

F{t) = [ /(f) dt' (Al) 

Jo 

the quantity to be determined. It is easy to prove that for U <t the exact identity 

F{t) = F(ti)e~^S + jf /(f) dt' (A2) 

holds, so that only the storage of values F(ti) and fiU) is necessary to evaluate F(ti+i). 
Note that equation (Al) is the solution of the differential equation 

f - /<«> - — ■ («) 

«t "Hag 

and this leads to an alternative way to compute F{t) by using (for instance) a finite difference 
integration scheme. In the integration of the hydrodynamical equations the evaluation of 
the quantity AF = F(t) dt over a time-step is also often required. From equations (Al) 
and (A2), simple algebra proves the exact relation 

AF = r Ug [Fit t ) - Fit l+1 )} + r lag f ' +1 f(t') dt'. (A4) 

Jti 

In the code, the function / is defined as the linear interpolation between the initial and final 
time over the time-step, so that the integrals (A2) and (A4) can be explicitly calculated with 
negligible computer time. 



5 Note that in CO01 the exponential factors in front of the integral in equation (B2) and inside the integral 
in equation (B3) are two typos. 
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Table V. Column 1 gives the adopted criterion to quantify fraction of the time spent in the 
high-luminosity states of accretion, in terms of the emitted SMBH bolometric luminosity and 
the current Eddington luminosity. In Column 2 we give the time percentage (calculated over 
the model bursting period) F Q n = A£ n/ Attesting, with At bursting = 5.5 Gyr (see Fig. [Q. 
Column 3 gives the time percentage of obscuration (more than 2 mag in the rest-frame UV) 
calulated from the ratio uv /0.2Lbh over the time interval AtoN- Finally, Fqso = -Fon x 
(1-Column 3) gives the fraction of the bursting period during which there is less than 2 mag 
of obscuration in the rest frame UV or less than 1.2 mag of obscuration in the rest-frame 
optical. This corresponds roughly to X-ray obscuration of ~ 10 22 atoms cm -2 . 



Lbr/ -^Edd 


-FoN 


At bsc/ AtoN 


Fqso 


> 0.100 


0.09% 


80% 


0.018% 


> 0.030 


1.51% 


94% 


0.096% 


> 0.010 


9.78% 


48% 


5.1% 


> 0.003 


22.48% 


36% 


14.4% 


> 0.001 


28.61% 


32% 


19.5% 
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